      Subroutine Aver(FIELD, GX)
      include "param.fi"
      include "common.fi"
      REAL FIELD(NTOT3T), GX(NTOT3T)
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance

      Do IB = 1, NTOT3T
        GX(IB) = FIELD(IB)
        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        I = HASHIND(II + DELI)
        Do DELJ = -1, 1
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        J = HASHIND(JJ + DELJ)
        Do DELK = -1, 1
        K = HASHIND(KK + DELK)
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC(X1(IB2),X2(IB2),X3(IB2),
     &                          X1(IB),X2(IB),X3(IB),
     &                          TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX + TMPY*TMPY + TMPZ*TMPZ)
            If (Distance.LE.HSUP .AND. Distance.GT.1E-10
     &         .AND. IB.NE.IB2) Then
              If (Shepard.EQ.0) Then
                GX(IB) = 0.5* (FIELD(IB2)-FIELD(IB)) * Weight(IB2)
     &                 * Kernel(Distance) + GX(IB)
              Else
                GX(IB) = 0.5* (FIELD(IB2)-FIELD(IB)) * Weight(IB2)
     &                 * Kernel(Distance) / SheCoef(IB) + GX(IB)
              Endif  
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
      Enddo

      End Subroutine
C==================================================================
      Subroutine Grad(FIELD, GX, GY, GZ)
      include "param.fi"
      include "common.fi"
      REAL FIELD(NTOT3T), GX(NTOT3T), GY(NTOT3T), GZ(NTOT3T)
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance
      REAL mtmp(2,2)

      Do IB = 1, NTOT3T
        GX(IB) = 0.0
        GY(IB) = 0.0
        GZ(IB) = 0.0

        If (tagCorr.EQ.CorrM1) Then
          Call GradCorrM1(IB, mtmp)
        Endif

        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        I = HASHIND(II + DELI)
        Do DELJ = -1, 1
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        J = HASHIND(JJ + DELJ)
        Do DELK = -1, 1
        K = HASHIND(KK + DELK)
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC(X1(IB2),X2(IB2),X3(IB2),
     &                          X1(IB),X2(IB),X3(IB),
     &                          TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX+ TMPY*TMPY+ TMPZ*TMPZ)
            If (Distance.LE.HSUP .AND. Distance .GT. 1E-10 
     &          .AND. IB.NE.IB2) Then
                Call KernelGradient(TMPX, TMPY, TMPZ, KX, KY, KZ)

              If (tagCorr.EQ.CorrM1) Call MatrixProd(mtmp, KX, KZ)

              If (Shepard.EQ.0) Then
                GX(IB) = GX(IB)-KX *(FIELD(IB2)-FIELD(IB)) * Weight(IB2)
                If (DIMM.EQ.DIMM3D)
     & GY(IB) = GY(IB)-KY *(FIELD(IB2)-FIELD(IB)) * Weight(IB2)
                GZ(IB) = GZ(IB)-KZ *(FIELD(IB2)-FIELD(IB)) * Weight(IB2)
              Else
                GX(IB) = GX(IB)-KX *(FIELD(IB2)-FIELD(IB)) * Weight(IB2)
     &                 / SheCoef(IB)
                If (DIMM.EQ.DIMM3D)
     & GY(IB) = GY(IB)-KY *(FIELD(IB2)-FIELD(IB)) * Weight(IB2)
     &                 / SheCoef(IB)
                GZ(IB) = GZ(IB)-KZ *(FIELD(IB2)-FIELD(IB)) * Weight(IB2)
     &                 / SheCoef(IB)
              Endif  
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
      Enddo

      End Subroutine
C==================================================================

      Subroutine Calc_ShepardCoef
      include "param.fi"
      include "common.fi"
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance

      Do IB = 1, NTOT3T
        SheCoef(IB) = 0.0
        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        I = HASHIND(II + DELI)
        Do DELJ = -1, 1
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        J = HASHIND(JJ + DELJ)
        Do DELK = -1, 1
        K = HASHIND(KK + DELK)
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC(X1(IB2),X2(IB2),X3(IB2),
     &                          X1(IB),X2(IB),X3(IB),
     &                          TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX+ TMPY*TMPY+ TMPZ*TMPZ)
            If (Distance.LE.HSUP) Then
                SheCoef(IB) = SheCoef(IB) +
     &                        Kernel(Distance) * Weight(IB2)
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
      Enddo

      End Subroutine

C==================================================================

      Subroutine Calc_dRho_dT
      include "param.fi"
      include "common.fi"
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance, mtmp(2,2)

      Do IB = 1, NTOT3T
        Temp1(IB) = 0.0
        If (tagCorr.EQ.CorrM1) Then
          Call GradCorrM1(IB, mtmp)
        Endif
        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        I = HASHIND(II + DELI)
        Do DELJ = -1, 1
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        J = HASHIND(JJ + DELJ)
        Do DELK = -1, 1
        K = HASHIND(KK + DELK)
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC(X1(IB2),X2(IB2),X3(IB2),
     &                          X1(IB),X2(IB),X3(IB),
     &                          TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX+ TMPY*TMPY+ TMPZ*TMPZ)
            If (Distance.LE.HSUP .AND. Distance .GT. 1E-10 
     &          .AND. IB.NE.IB2) Then
                Call KernelGradient(TMPX, TMPY, TMPZ, KX, KY, KZ)
              If (tagCorr.EQ.CorrM1) Call MatrixProd(mtmp, KX, KZ)
              If (Shepard.EQ.0) Then
                Temp1(IB) = Temp1(IB) -
     &            KX * (U1(IB2)-U1(IB)) * Weight(IB2) * Rho(IB2)
                If (DIMM.EQ.DIMM3D)
     &           Temp1(IB) = Temp1(IB) -
     &            KY * (U2(IB2)-U2(IB)) * Weight(IB2) * Rho(IB2)
                Temp1(IB) = Temp1(IB) -
     &            KZ * (U3(IB2)-U3(IB)) * Weight(IB2) * Rho(IB2)
              Else
                Temp1(IB) = Temp1(IB) -
     &            KX * (U1(IB2)-U1(IB)) * Weight(IB2) * Rho(IB2)
     &            / SheCoef(IB)
                If (DIMM.EQ.DIMM3D)
     &           Temp1(IB) = Temp1(IB) -
     &            KY * (U2(IB2)-U2(IB)) * Weight(IB2) * Rho(IB2)
     &            / SheCoef(IB)
                Temp1(IB) = Temp1(IB) -
     &            KZ * (U3(IB2)-U3(IB)) * Weight(IB2) * Rho(IB2)
     &            / SheCoef(IB)
              Endif
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
      Enddo

      End Subroutine
C==================================================================

      Subroutine Calc_dU_dT
      include "param.fi"
      include "common.fi"
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance, mtmp(2,2)

      Do IB = 1, NTOT3T
        Temp1(IB) = 0.0
        Temp2(IB) = 0.0
        Temp3(IB) = 0.0
        If (tagCorr.EQ.CorrM1) Then
          Call GradCorrM1(IB, mtmp)
        Endif
        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        I = HASHIND(II + DELI)
        Do DELJ = -1, 1
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        J = HASHIND(JJ + DELJ)
        Do DELK = -1, 1
        K = HASHIND(KK + DELK)
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC(X1(IB2),X2(IB2),X3(IB2),
     &                          X1(IB),X2(IB),X3(IB),
     &                          TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX+ TMPY*TMPY+ TMPZ*TMPZ)
            If (Distance.LE.HSUP .AND. Distance .GT. 1E-10 
     &          .AND. IB.NE.IB2) Then
                Call KernelGradient(TMPX, TMPY, TMPZ, KX, KY, KZ)
!                If (IB.EQ.1) Write(*,*) 
!     &            P(IB)/Rho(IB)/Rho(IB) , P(IB2)/Rho(IB2)/Rho(IB2)
!     &            ,Weight(IB2), KX
              If (tagCorr.EQ.CorrM1) Call MatrixProd(mtmp, KX, KZ)
              If (Shepard.EQ.0) Then
                Temp1(IB) = Temp1(IB) -
     &            (-P(IB)/Rho(IB)/Rho(IB) + P(IB2)/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KX
                If (DIMM.EQ.DIMM3D)
     &          Temp2(IB) = Temp2(IB) -
     &            (-P(IB)/Rho(IB)/Rho(IB) + P(IB2)/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KY
                Temp3(IB) = Temp3(IB) -
     &            (-P(IB)/Rho(IB)/Rho(IB) + P(IB2)/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KZ
              Else
                Temp1(IB) = Temp1(IB) -
     &            (-P(IB)/Rho(IB)/Rho(IB) + P(IB2)/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KX
     &            / SheCoef(IB)
                If (DIMM.EQ.DIMM3D)
     &          Temp2(IB) = Temp2(IB) -
     &            (-P(IB)/Rho(IB)/Rho(IB) + P(IB2)/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KY
     &            / SheCoef(IB)
                Temp3(IB) = Temp3(IB) -
     &            (-P(IB)/Rho(IB)/Rho(IB) + P(IB2)/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KZ
     &            / SheCoef(IB)
              Endif
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
      Enddo

      End Subroutine

C==================================================================
      Subroutine Calc_VisTerm
      include "param.fi"
      include "common.fi"
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance

      Do IB = 1, NTOT3T
        Temp1(IB) = 0.0
        Temp2(IB) = 0.0
        Temp3(IB) = 0.0
        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        I = HASHIND(II + DELI)
        Do DELJ = -1, 1
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        J = HASHIND(JJ + DELJ)
        Do DELK = -1, 1
        K = HASHIND(KK + DELK)
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC(X1(IB2),X2(IB2),X3(IB2),
     &                          X1(IB),X2(IB),X3(IB),
     &                          TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX+ TMPY*TMPY+ TMPZ*TMPZ)
            If (Distance.LE.HSUP) Then
                KX = KernelLaplace(Distance)
              If (Shepard.EQ.0) Then
                Temp1(IB) = Temp1(IB) +
     &            KX *(U1(IB2)-U1(IB)) * Weight(IB2) * NU
                If (DIMM.EQ.DIMM3D)
     &           Temp2(IB) = Temp2(IB) +
     &            KX *(U2(IB2)-U2(IB)) * Weight(IB2) * NU
                Temp3(IB) = Temp3(IB) +
     &            KX *(U3(IB2)-U3(IB)) * Weight(IB2) * NU
              Else
                Temp1(IB) = Temp1(IB) +
     &            KX *(U1(IB2)-U1(IB)) * Weight(IB2) * NU
     &            / SheCoef(IB)
                If (DIMM.EQ.DIMM3D)
     &           Temp2(IB) = Temp2(IB) +
     &            KX *(U2(IB2)-U2(IB)) * Weight(IB2) * NU
     &            / SheCoef(IB)
                Temp3(IB) = Temp3(IB) +
     &            KX *(U3(IB2)-U3(IB)) * Weight(IB2) * NU
     &            / SheCoef(IB)
              Endif
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
      Enddo

      End Subroutine

C==================================================================
      Subroutine Calc_NLTerm
      include "param.fi"
      include "common.fi"

      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance

      If (DIMM.EQ.DIMM3D) Then
        Do I = 1, NTOT3T
          Temp1(I) = - U1(I) * DU1DX(I) 
     &             - U2(I) * DU1DY(I) - U3(I) *DU1DZ(I)
          Temp2(I) = - U1(I) * DU2DX(I) 
     &             - U2(I) * DU2DY(I) - U3(I) *DU2DZ(I)
          Temp3(I) = - U1(I) * DU3DX(I) 
     &             - U2(I) * DU3DY(I) - U3(I) *DU3DZ(I)
        Enddo
      Else
        Do I = 1, NTOT3T
          Temp1(I) = - U1(I) * DU1DX(I) 
     &             - U3(I) *DU1DZ(I)
          Temp3(I) = - U1(I) * DU3DX(I) 
     &             - U3(I) *DU3DZ(I)
        Enddo
      Endif

      End Subroutine

